setwd("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207")
filea=read.csv("all.FDR.sig.at.least.one.add.direction.same.diff.csv",header=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]

file=read.csv("at.least.one.AShM.in.DC.add.BF.beta0.add.CCHC.csv",header=T)
file$id=paste(file$Chr,file$Start,sep=":")
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]
eqtl=fread("/mnt/md1200/5/zhaocunyou/wangzhongju/eQTL_analysis_by_wangzhongju/eQTL_data/GTEx_v6p_eQTL_for_lqy.txt",header=T,sep="\t")
eqtl1=tidyr::separate(eqtl,SNP,into=c("Chr","Pos","Alt","Ref","Other"),sep="_")
eqtl1$Chr=paste0("chr",eqtl1$Chr)
eqtl1$unitID=paste(eqtl1$Chr,eqtl1$Pos,eqtl1$Alt,eqtl1$Ref,sep=":")
rt=merge(file,eqtl1,by="unitID",all.x=T)
rt1=rt[rt$pattern.not.rm.dupl.num.DC>1,]
rt2=rt1[rt1$BF_in_DC>1,]
rt3=rt1[rt1$BF_in_DC>10,]

rt3=rt3[!is.na(rt3$P),]
length(unique(rt3$unitID))

tidyr::unit
library(data.table)
gwas.scz1=fread("/mnt/md1200/5/zhaocunyou/wangzhongju/GWAS_data_by_wangzhongju/GWAS/ckqny.scz2snpres",header=T,sep="\t")#精分的GWAS data
gwas.scz1=gwas.scz1[gwas.scz1$p<0.05,]
gwas.scz1$id=paste(gwas.scz1$hg19chrc,gwas.scz1$bp,sep=":")
gwas.scz1=data.frame(gwas.scz1$id,group.gwas=rep("scz1",dim(gwas.scz1)[1]))

gwas.scz2=fread("/mnt/md1200/5/zhaocunyou/wangzhongju/GWAS_data_by_wangzhongju/GWAS/clozuk_pgc2.meta.sumstats.ma",header=T,sep="\t")#精分的GWAS data
gwas.scz2=gwas.scz2[gwas.scz2$p<0.05,]
gwas.scz2$chr=paste("chr",gwas.scz2$chr,sep="")
gwas.scz2$id=paste(gwas.scz2$chr,gwas.scz2$pos,sep=":")
gwas.scz2=data.frame(gwas.scz2$id,group.gwas=rep("scz2",dim(gwas.scz2)[1]))

gwas.bd=fread("/mnt/md1200/5/zhaocunyou/wangzhongju/GWAS_data_by_wangzhongju/GWAS/daner_PGC_BIP32b_mds7a_0416a",header=T,sep="\t")#双相的GWAS data
gwas.bd=gwas.bd[gwas.bd$P<0.05,]
gwas.bd$CHR=paste("chr",gwas.bd$CHR,sep="")
gwas.bd$id=paste(gwas.bd$CHR,gwas.bd$BP,sep=":")
gwas.bd=data.frame(gwas.bd$id,group.gwas=rep("bd",dim(gwas.bd)[1]))

names(gwas.scz1)=c("id","group.gwas.scz1")
names(gwas.scz2)=c("id","group.gwas.scz2")
names(gwas.bd)=c("id","group.gwas.bd")

rt3=merge(rt3,gwas.scz1,by="id",all.x=T)
rt3=merge(rt3,gwas.scz2,by="id",all.x=T)
rt3=merge(rt3,gwas.bd,by="id",all.x=T)




rt3$eqtl.no.na.num=rowSums(!is.na(rt3[,178:179]),na.rm=T)/2
rt3$gwas.no.na.num=rowSums(!is.na(rt3[,184:186]),na.rm=T)
write.csv(rt3,"200ASH.add.eQTL.GWAS.csv",quote=F,row.names=F)

###20201223做GWAS的富集分析。case是807 psyASH，control是117K-807,在服务器上操作
cd /mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201222

names(gwas.scz1)=c("id","group")
names(gwas.scz2)=c("id","group")
names(gwas.bd)=c("id","group")
gwas=rbind(gwas.scz1,gwas.scz2,gwas.bd)
gwasid=unique(as.character(gwas$id))

con1=filea[!filea$id %in% file2$id,]
c=length(intersect(as.character(con1$id),as.character(gwasid)))
d=dim(con1)[1]-c

result=data.frame(matrix(,nrow = 3,ncol = 4))
names(result)=c("Term","overlap.num","OR","P.value")
result[1,]=c("control.all.ASH.rm.psyASH",c,NA,NA)

a=length(intersect(as.character(file2$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,807-a,d),nrow = 2,byrow = T))
result[2,]=c("807",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file3$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,200-a,d),nrow = 2,byrow = T))
result[3,]=c("200",a,tmp$estimate,tmp$p.value)
#结果不显著，因此换case
> result
                       Term overlap.num               OR           P.value
1 control.all.ASH.rm.psyASH       24323             <NA>              <NA>
2                       807         190 1.16330640394036 0.075016400387453
3                       200          46 1.12839198224925  0.48634572461143


###
allsnp=fread("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/all.heter.genotype.snp",header=T,sep="\t")

con=allsnp[!allsnp$id %in% as.character(filea$id),]		###remove the ASH
conid=as.character(con$id)
conid.random=sample(conid,150000,replace = FALSE)
con1=con[con$id %in% conid.random,]
c=length(intersect(as.character(con1$id),as.character(gwasid)))
d=dim(con1)[1]-c

result=data.frame(matrix(,nrow = 7,ncol = 4))
names(result)=c("Term","overlap.num","OR","P.value")
result[1,]=c("control.all.heter.rm.ASH",c,NA,NA)

a=length(intersect(as.character(filea$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,117012-a,d),nrow = 2,byrow = T))
result[2,]=c("117K",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(filea1$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,61725-a,d),nrow = 2,byrow = T))
result[3,]=c("61K",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,53425-a,d),nrow = 2,byrow = T))
result[4,]=c("53425",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file1$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,8544-a,d),nrow = 2,byrow = T))
result[5,]=c("8544",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file2$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,807-a,d),nrow = 2,byrow = T))
result[6,]=c("807",a,tmp$estimate,tmp$p.value)

a=length(intersect(as.character(file3$id),as.character(gwasid)))
tmp=fisher.test(matrix(c(a,c,200-a,d),nrow = 2,byrow = T))
result[7,]=c("200",a,tmp$estimate,tmp$p.value)

write.csv(result,"./GWAS.random.con.1.statis.csv",quote=F,row.names = F)
###算meta.p
library(meta)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201223.GWAS")
sel=list.files(pattern = "GWAS.random.con")
f1=read.csv(sel[1],header=T)
f2=read.csv(sel[2],header=T)
f3=read.csv(sel[3],header=T)
info.data=data.frame(Term=c("117K","61K","53425","8544","807","200"),
                     unique.num=c(117012,61725,53425,8544,807,200))

result=data.frame(matrix(,nrow = 6,ncol = 5))
names(result)=c("Term","OR","P.value","upper","lower")
for(i in 1:6){
case.overlap.num=as.numeric(rep(f1[i+1,"overlap.num"],3))
case.total.num=as.numeric(rep(info.data[i,"unique.num"],3))
con.overlap.num=as.numeric(c(f1[1,"overlap.num"],f2[1,"overlap.num"],f3[1,"overlap.num"]))
con.total.num=as.numeric(rep(150000,3))

tdata=data.frame(case.overlap.num,case.total.num,con.overlap.num,con.total.num)
metaor3<-metabin(case.overlap.num,case.total.num,con.overlap.num,con.total.num,data=tdata,sm="OR")
result[i,1]=info.data[i,1]
result[i,2]=OR=exp(metaor3$TE.fixed)
result[i,3]=metaor3$pval.fixed
result[i,4]=upper=exp(metaor3$upper.fixed)
result[i,5]=lower=exp(metaor3$lower.fixed)
}
result$FDR=p.adjust(result$P.value,method = "BH")
write.csv(result,"./meta.p.GWAS.statis.csv",quote=F,row.names = F)
